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ABSTRACT 

We consider the problem of automatically (and robustly) isolating and extracting 
information about waves and oscillations observed in EUV image sequences of the 
solar corona with a view to near real-time application to data from the Atmospheric 
Imaging Array (AIA) on the Solar Dynamics Observatory (SDO). We find that a sim- 
ple coherence / travel-time based approach detects and provides a wealth of informa- 
tion on transverse and longitudinal wave phenomena in the test sequences provided by 
the Transition Region and Coronal Explorer (TRACE). The results of the search are 
pruned (based on diagnostic errors) to minimize false-detections such that the remain- 
der provides robust measurements of waves in the solar corona, with the calculated 
propagation speed allowing automated distinction between various wave modes. In 
this paper we discuss the technique, present results on the TRACE test sequences, and 
describe how our method can be used to automatically process the enormous flow of 
data (^lTb/day) that will be provided by SDO/AIA after launch in late 2008. 



Since the launch of the Transition Region and Coronal Explorer (TRACE; lHandy et a/.l ll999) 
the community has invested a great deal of effort on the identification and analysis of longitudi- 
nal and transverse wave phenom ena in loop structures seen in its EUV images of the corona (see 



Nakariakov and Verwichte 20051 for a good overview). The great interest in finding and character- 



izing coronal waves is driven by the promise of coronal seismology which could lead to the deter- 
mination of otherwise inaccessible physical properties of th e solar atmosphere by studying p hase 



speeds, amplitudes, dissipation, etc. of the observed waves (INakariakov and Verwichtel l2005). 



There are many unresolved issues in coronal seismology: for example, it is unclear why 
only a subset of coronal loops show transverse oscillations in the wake of a flare or CME {e.g., 
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Aschwanden et a/.ll2002l). which physical processes dominate the damping of these transverse os 



dilati ons {e.g. , Ofman and Aschwandenll2002 ; Arregui et a/]|2007 ; Terradas. Andries. and Goossens 



2007 ) . why only a subset of coronal loops show propagating slow-mode oscillations (e.g. Jde Moortel. Ireland, and 
200o]; Ide Moortel and Rosner!l2~007l) . how the propagation of slow-mode waves at different temper- 



atures can be used to probe the thermal structure and loop scale heights (King et al. 



2003), or how 



lower-atmospheric oscillat i ons generally leak into the corona (e.g.. iDe Pontieu. Erdelyi. and De Moortel 
20051 : Ijefferies et a/lbood : iKhomenko efaZlbood) ? 



Many of these issues require more extensive study with much larger statistical samples. It 
is surprising that the number of wave and oscillation detections in the above papers is only of or- 
der several dozen for both longitudinal and transverse oscillations. This is a very small number 
compared to the volume of TRACE data given that flares occur often and "wave-leakage" from 
the solar interior to the lower portions of the outer solar atmosphere is apparently quite abundant 
(be Pontieu. Erdelvi. and De Moortelll2005l : ljefferies et a/.lbood : IKhomenko et a/.ll2008h . The dis- 
crepancy between the plethora of observed wave driven phenomena in the chromosphere and tran- 
sition region, and the relative paucity of waves seen in the corona partly motivates the effort dis- 
cussed herein. Is the small number of coronal wave observations caused by the fact that all of the 
TRACE observations of waves have been found as a result of a "manual" search, i.e., limited by 
human patience and detection skills? Or are there actually very few locations in the corona where 
waves can be observed? 

One of our aims is to develop an automated wave-detection algorithm to help expand the 
number of known oscillations, and build up significant statistics on how frequent coronal oscilla- 
tions occur. Such larger statistical samples are crucial to address many of the unresolved issues 
in coronal seismology, and will be within reach with the Atmospheric Imaging Assembly (AIA) 
after launch of the Solar Dynamics Observatory (SDO) by 2009. This instrument will provide 
full-disk coronal imaging at significantly increased signal-to-noise in eight coronal wavelengths at 
ten-second cadence continuously, producing an enormous data flow of ^lTb day -1 . The potential 
for significantly increased detection of coronal oscillations and waves with such a data rate and 
quality are very high. However, AIA's enormous data rate renders business-as-usual approaches 
that involve manually looking through data sets and individually flagging events unfeasible. Since 
it is unpractical to transfer the full volume of AIA data to remote users, it is critical that auto- 
mated algorithms be developed that can search for wave-like phenomena and automatically flag 
temporal and spatial regions of interest for later, more detailed analysis. There have been a few 
papers recently that have started addressing some of the issues that need to be resolved to enable 
full optimal use of the AIA data, such a s automated detection of locations with significant oscilla- 
tory po wer (Nakariakov and Kingll2007|) . or semi- automated detection of propagating and standing 
waves (|Sych and Nakariakovll2008|) . The approach presented in this paper may be the first that can 
reliably and automatically distinguish between longitudinal and transverse oscillations at the same 



-3- 



time as rejecting most false positive detections of oscillations. Our approach is base d on the tech 



niqu e used to analyze and detect coronal Alfven waves with the CoMP instrument (|Tomczyk et al 



20071) 



In the following section we discuss the method developed, illustrating its application on ex- 
ample TRACE datasets (see Sect. |2]) which are known to show a host of transverse and longitudinal 
wave phenomena and are therefore good for testing our approach. In Sect. [3] we explore the results 
of the analysis on the test data and discuss the method employed to minimize false detections. In 
Sect. ID we compare some of our results with previous analyses of the TRACE datasets we have 
used, and discuss how our method could be used in a practical manner for SDO/AIA data. 



1. Method 



We adapt the phase travel-time analysis (e.g. Jjefferies et a/.lll 99411 1 9971 : LFinsterle et a/.ll2004 : 
Mcintosh. Fleck, and Tarbellll2004|) to isolate the propagation characteristics of the wave modes 
observed in the TRACE image sequences time-series as a function of frequency. This technique is 
performed as a Fourier comparison of time-series in neighboring image pixels to extract informa- 
tion about the cross-correlation, coherence, cro s s-po wer and phase dependence of the neighbor- 
ing signals. As demonstrated in iTomczyk et all (|2007|) this is a particularly simple, but powerful 
technique to track coronal waves, albeit in the plane-of-the-sky. The technique has the inherent 
flexibility to be employed rapidly in the spectral domain producing robust results. 

The first step in the travel-time analysis involves computing the multi dimensional Fourier 
transform of the EUV image cube (I(x, y, t)) in the third dimension which gives a complex cube 
If(x, y, f) that now has frequency / [units: mHz] as its third dimension. Before calculating the 
Fourier transform we apply a Planning window on the timeseries to reduce the impact of temporal 
trends in the short data sequences. 

The analysis involves comparisons of the intensity timeseries of our pixel of interest (x, y) 
and that of a pixel within a square region (of dimensions dx x dy) around the central pixel. At 
this point the analysis diverges into t he spectral (Sect. II. II) an d temp oral domains (Sect. 11.21) where 
the former has been developed by ITomczyk. and Mcintosh! (|2008l) to supersede the analysis of 
Tomczyk et al.\ (|2007|) . We have chosen the range to be ±10 pixels in both spatial dimensions. 
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1.1. Spectral Domain Analysis 

For the range of pixels [x ± dx : y ± dy] we compute the cross-spectrum 

cs RB (f) = i* R (f)*r B (fy (i) 

between the Fourier transformed timeseries at the reference pixel I R (f) for any other pixel in the 
analysis box I B (f), with I(fY being the complex conjugate of /(/). We then smooth the cross- 
spectrum with a five-point boxcar smoothing function in frequency to reduce the contribution of 
noise in the computation. In the spectral domain we wish to rapidly compute wave diagnostics in 
several characteristic frequency filters simultaneously and for the sake of simplicity we will use 
Gaussian filters of the form G(f , 5f) where f is the central frequency and 8f is the width of the 
Gaussian filter chosen which, for this paper, is /o/10. So, for each filter selected we can simul- 
taneously compute the weighted signal cross-power [WCP RB (f )], phase [4>R B (fo)] , coherence 
[C RB (fo)] and phase travel-time [T RB (f )] for the pixel of interest that are given by the following 
expressions: 

WCPrbUo) = G(f ,5f)\\CS RB (f)\\ (2) 



>RB 



(Jo) 



G(/ °'' /)tan {Re{C RB (f)}) (3) 

CM) = G(f , 8f) (cs^0c%m) (4) 

-Irb{Jo) = 7 — y?) 

Repeating this step for all of the neighboring pixels results in the construction of lOx 10 pixel maps 
of the three quantities for each filter frequency, / . 

The panels in Figured] show coherence, cross-power, phase and phase travel-time maps from 
the spectral analysis of one pixel (coordinates 181, 97) in the TRACE 171 A dataset from 1 July 
1998 12:03, - 1 July 1998 13:02 UTC (see below) at a filter frequency of /„ = 5.5mHz. In particular , 



we see that the coherence map (panel A), like that shown in the approach of iTomczyk et al.\ (|2007|) . 
has a region, or "island", where the spectral coherence is large (>0.6) at this frequency across 
numerous neighboring pixels. It is the properties of the strong coherence island and its mappings in 
phase-time and cross-power that form the cornerstone of the analysis presented. They can provide a 
quick and reliable estimate of the propagation direction and phase-speed of the disturbance present, 
from which reasonable errors can be determined. 

After isolating the relevant high coherence pixels in the phase-time diagnostic map, we com- 
pute a distance-minimizing linear fi{]to the cluster of points in the island and calculate a distance d 



'in this case we seek the best fitting straight line through the cluster of points in the island that also minimizes the 
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along the resulting straight line as well as an angle (and the error on the angle) at which the distur- 
bance is moving (relative to North - South in the coherence map). We also use the cluster of pixels 
to extract a measure of two coherent length- scales, one that is parallel to the measured propagating 
direction (the coherence length), and one that is perpendicular to it (the coherence width). These 
last diagnostics may contain physical information about the process exciting the disturbance and 
so we retain them at little computational cost. 

We then compute the phase-speed of the disturbance and its associated error by forming a 
scatter plot of distance d (from the reference pixel) versus measured phase-time (see, e.g., Fig- 
ure [2]). The gradient of the least-squares linear fit to these data provides an estimate of the phase 
speed of the disturbance, errors in the gradient are determined by assessi ng a range of possible val 



ues th at can achieve a similar quality of fit to the data (see, e.g. Sect. 6 of iBevington and Robinson 



20031). We repeat the previous step for each filter required unless there are fewer than 6 pixels in 



the island. 

This process is then repeated for all of the other pixels in the original image with the result 
that we develop images of the apparent phase- speed, angle, (their respective errors) and additional 
measures of the island-mean coherence and cross-power (see, e.g, Figure[5]and Sect. [3]). 



1.2. Temporal Domain Analysis 



At the cost of a considerable increase in computation time, we can substitute the phase-time 
and phase-speed measurements by performing their calculation in the temporal domain. This step 
involves the cross-correlation of the Fourier filtered timeseries for all of the pixels in the high 
coherence island relative to the reference pixel. The cross-correlation function at each neighboring 
pixel is a Gabor wavelet that, when fitted , yields information about the phase (and energy carrying 
group) travel-times of the disturbance (IFinsterle et all 120041) . In practice we cannot execute the 
temporal domain analysis for SDO/AIA as the computational cost is prohibitive for a method that 
is essentially mathematically equivalent to the one presented above. However, we do note that 
the process is illustrative of the general technique, the quantities and measures developed for the 
analysis and so we have chosen to provide an example. 



distance between the points and that straight line (see, e.g. Bevington and Robinson! 2003 ). We have found that this 
approach is not subject to singularities for clusters of points that are orthogonal, or nearly orthogonal, to the horizontal 
as is the case in standard least-squares fits when the dependent and independent variables are interchangeable or both 
have measurement errors. 
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An example of two-stepj cross-correlation fitting is shown in Figure [3l First, we fit the 
envelope of the function to extract group information, after which we fit the peak of the cross- 
correlation nearest to the lag origin to obtain phase information. The solid black line shows the 
cross-correlation of the reference timeseries used to create Figure Q] and that from another pixel 
(coordinate [181, 91] in the figure) as a function of lag in seconds. The group behavior is deter- 
mined by fitting a Gaussian (dot-dashed line) to the envelope of the cross-correlation function that 
is outlined by the absolute values of the function maxima (red triangles) - the center of the fitted 
Gaussian measures the group travel-time (or energy transport time). The phase travel-time can then 
be extracted by fitting the peak of the cross-correlation function nearest to the origin. Using either 
a parabolic of Gaussian fit we can extract phase-times (the center of the Gaussian or parabolic 
fit) down to an accuracy of about one tenth of the image cadence. For the case considered the 
vertical dashed lines in panel C of Figure [3] show the positions of the estimated group travel-time 
(red; -4.17 minutes) and phase travel-time (green; -0.527 minutes. For good measure, we note that 
the latter value is consistent with the values in Figures Q] and demonstrates that the disturbance is 
propagating left to right which is also consistent with the movie of Figure HI) respectively. 

Repeating the steps used in the spectral analysis to compute the phase-speed we perform a 
least- squares linear- fit to the scatter plot of pixel distance and measured phase-travel times. We 
note that it is also possible to investigate the group-speed of the disturbance - by simply replacing 
the phase travel-times with the measured group travel-times and repeating the fit to the scatter. 



1.3. Improving Algorithm Efficiency 

The algorithm described in Sects. [TTTl and |1.2| is quite slow, since it calculates cross-correlations 
or cross-spectral properties for 100 neighboring points at each location of the datacube. We have 
improved the efficiency of the algorithm by a large factor (typically of order 20) in the following 
way. For each location, we define a map of locations for which cross-correlations with the central 
pixel [xo, yo] need to be calculated. At the first iteration this map contains only the 4 immediate 
neighbors of the central pixel. We then determine which of these locations has a coherence larger 
than the cutoff (0.6 in our case). If none of these locations has a coherence larger then the cutoff, 
the algorithm ends the calculations for pixel [xq, t/oL an d moves on to [x + 1, y ]. If one or more 
locations shows a coherence larger than the cutoff, we repeat the iteration by calculating a new 
map of points for which cross-correlations need to be calculated. This map will now contain the 



2 Typically we would fit both phase and group information simultaneously, as is done in lFinsterle et al\ d2004l) and 
the other references provided above, but for the data at hand, where often the group envelope is not well defined, we 
fit them separately. 
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immediate neighbors of the high coherence locations calculated in the previous round. We exclude 
from this map all locations that have already been calculated. After this the algorithm again checks 
whether any of the new points has a coherence larger than the cutoff, continuing until no new 
locations are above the cutoff. In this fashion we can usually reduce the number of calculations 
from 100 to around 5, since many locations show little coherence with their neighbors, and only 
a select few locations show significant coherence with their neighbors. This increase in efficiency 
has a critical impact on the usefulness of the algorithm for searching through large datasets, such 
as those that are expected from SDO/AIA. 



2. Test Data 

To illustrate the algorithm we have applied it on two different TRACE datasets. Dataset I was 
taken on 1 July 1998 from 12:03 to 13:02 UTC in 17lA, with a cadence of 31.5s. Dataset II was 
taken on 14 July 1998 from 12:45 to 13:42 UTC in 17lA during the so-called Bastille day flare 
and has a cadence of 72.3 seconds. The first dataset contains an active region with some slow- 
mode running waves in one of the loops associated with a sunspot. T he second dataset has been 



Aschwanden et a/.l ll999. 



2002) and contains the "classical 



the subject of several papers (see, e.g. 
example" of transverse oscillations. Before applying the algorithm, we performed the calibration 
steps included in trace_prep.pro (as part of the TRACE tree of solarsoft IDL), which includes the 
subtraction of dark current/pedestal, correction for the ccd gain and lumogen degradation, cor- 
rection of bad pixels and normalization of all images to remove differences in exposure time. In 
addition, we performed several iterations of cosmic ray removal using the trace_unspike.pro rou- 
tine which removes excessively bright spikes by comparing the data with previous and following 
timesteps and replacing the spike with a temporal and spatial average of the surrounding pixels. 
As a final step we derotate the data by using poly_2d.pro to perform sub-pixels shifts that are cal- 
culated using difLrot.pro. We also perform a 2x2 rebinning to reduce the TRACE data to a spatial 
resolution of 1 arcsecond. We have tested the impact of the rebinning extensively and found that 
rebinning to the spatial resolution does not lead to any loss of useful information about the os- 
cillations, and in fact improves the signal to noise of the sometimes weak oscillatory signals. Of 
course, rebinning the data naturally speeds up the algorithm by a factor of 4 compared to running 
it on the data at its original resolution. 

Figure |4] show snapshots of the raw and Fourier filtered timeseries of the 1 July 1998 (left) 
and July 14 1998 (right) datasets respectively. In each case, the upper left panel shows a snapshot 
of the intensity timeseries while other panels show the filtered timeseries with a Gaussian filter 
centered around 1.5 mHz, 3.5 mHz and 5.5 mHz respectively. In the left panel we see that the 
active region contains several sets of loops, with the right-most fan associated with a sunspot 
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being highlighted in a box. The loops in the southern part of the active region are associated 
with plage regions. We see that the sunspot fan contains a clear 180s signal which is visible as a 
dark spot in the 5.5 mHz snapshot (lower right) and an oscillatory signal in the same panel of the 
pr ovided movie. These kind s of oscillations have been described extensively in a series of papers 



by Ide Moortel et all (|2002aUb|) . In addition, the moss regions towards the north and middle of the 



active region sh ow a predominance of 1.5 and 3.5 mHz power. This is not surprising and has been 



noted before by lDe Pontieu. Tarbell. and Erdelyil (|2003|) who showed that the periodic obscuration 



of upper transition region emission by chromospheric jets plays a dominant role in the EUV (semi- 
periodicities. A very intriguing finding is that there is significant power in the 1.5 mHz filtered 
movies and snapshot in the plage-related loops towards the south of the active region. We will 
return to these low frequency oscillations in Sects. [3] and HI 

As a cautionary note we see that generally, these filtered movies show not only periodic sig- 
nals, but any signal that shows power at the chosen frequency. Because of the nature of the Fourier 
transform, this includes signals that are caused by abrupt changes in the brightness, such as cosmic 
rays or sudden loop brightenings (e.g., in the case of emerging flux reconnecting with overlying 
field, or in post-flare loop arcades). To exclude some of those signals when we run our algorithm 
on the data, we take care to use a Hanning window (reducing the influence of trends), and apply a 
mask that excludes from our calculations all pixel locations whose raw time series show deviations 
that are more than 3.5 times the standard deviation of the time series. This practically excludes 
many of the locations with cosmic rays and sudden brightenings. 

This approach is especially useful in the case of a flare where many locations show significant 
and sudden brightenings, such as is the case in dataset II, and shown in the right panel. The filtered 
movies and snapshots show many locations with significant oscillatory power, which is borne out 
by the unfiltered time series. The flare is associated with a coronal dimming, which leads to 
significant "wiggling" of the coronal loop structures, especially on the 5-10 minute time scales. 
Very prominent in these snapshots are also the post-flare loops in the middle of the active region. 
Most interesting is that the transversely oscillating loops show up very clearly in the movies and 
snapshots, especially at 3.5 mHz (lower left) and to a lesser degree at 5.5 mHz (lower right). 



3. Results 

3.1. Interpretation of Raw Results 

We present results of the spectral algorithm for both datasets. As we have stated above the 
results for the temporal code are very similar but run at a much lower speed so that, as such, it is not 
adequate for our ultimate goal of applying this technique to SDO/AIA. In principle our algorithm 
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produces for each location, and for each frequency we choose to filter on, the following average 
properties of the island of high coherence around each location: cross-spectral power, average 
coherence, phase speed of propagating signal, error on the phase speed, the angle of propagation 
of the disturbance, the correlation length and the correlation width. To reduce the amount of noise 
in the resulting maps of the active region, we reject locations for which the island of high coherence 
consists of fewer than 6 neighbors. 

The resulting maps of the properties described above show only locations where the signal is 
significantly correlated with at least 5 neighbors (Figures [51 [6] and [7] for the frequencies 1.5, 3.5 
and 5.5 mHz for dataset I, and Figures [8l [9]and[10]for dataset II). What is striking in these maps 
is the large amount of locations that have (oscillatory) correlation with their neighbors. However, 
since our algorithm is so sensitive, many of these locations actually show "false positives" in which 
a very small (mostly instrumental or spurious) oscillatory signal is correlated with a few neighbor- 
ing pixels. Nevertheless, some interesting physical results immediately show up that confirm the 
findings of our cursory analysis of the Fourier- filtered time series in Section 3. For example, the 
3.5 mHz map of phase speeds of dataset I (Figure [6t is dominated by locations where mo ss oc- 



curs, as can be expected from previous analyses (e.g.. De Pontieu. Tarbell. and ErdeTyill2003h . The 
phase speeds are quite variable and of order 10-50 km s _1 in the moss regions, which is perhaps 
not surprising since the correlations with neighboring pixe ls are most likely caused by chromo 



spheric jets (with sim ilar velocities to those we find; (e.g. JPe Pontieu. Tarbell. and Erdelyill2003 



Hansteen et al.\ 12006) moving across TRACE pixels. The 3.5 mHz map also shows evidence of 
oscillatory power in the bottom of some coronal loops (e.g., around 1 10, 90) with phase speeds of 
order 50- 100 km s 1 . 

The 5.5 mHz map of dataset I (Figure [7]) is much more sparsely populated with only a few 
major concentrations of significant oscillatory power towards the lower right of the map. These 
locations are exactly where we find (through visual examination) oscillations in the original and 
Fourier- filtered time series: our algorithm nicely locates the sunspot oscillations at 180 second 
periods (around 180, 90). The phase speeds here are of order 5 0-150 km s" 1 , which is within the 



range that is expected for slow-mode propagating disturbances (|de Moortel et a/.ll200 2a.b) 



The 1.5 mHz results (Figure [5]) show large parts of the southern coronal loops with significant 
oscillatory disturbances that propagate along the loops with phase speeds that are similar and of 
order 50 - 150 km s -1 . If these oscillations are real (see Sect. [5]), they are interesting, since they are 
of much longer periods than have been previously reported (typically five minutes in plage-related 
loops). We have verified that the signals are not introduced by slow pointing drifts on timescales 
of ten minutes by performing the same analysis using data that is co-aligned with various co- 
alignment strategies (see Sect. |5]). In addition, co-alignment issues would lead to similar signals in 
similarly oriented loops, e.g., the loops at the northern edge of the sunspot fan. Our algorithm does 
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not find any sign of significant propagating and oscillatory disturbances at 1.5 mHz in those loops. 

The results of dataset II are similarly interesting. First of all, the phase speeds shown in all 
frequencies are systematically higher in this dataset, with most being of order 150 - 500 km s _1 . 
This is compatible with the fact that this dataset is dominated by a flare, coronal dimming and the 
associated (mostly fast magneto-acoustic mode) waves and/or standing waves (transverse oscilla- 
tions) that have been well docu mented in previous work focusing on the Bastille day flare (e.g., 



Aschwanden et a/.N 19991 12002|) . The 3.5 mHz maps (Figure |U) of phase speed are dominated by 
locations of significant oscillatory power at the same locations where previous analysis (based on 
visual inspection) had found evidence of transversely oscillating loops (e.g., around 170, 200 or 
180, 70 or 240, 1 10). In addition, our algorithm finds a range of locations with oscillatory power 
that have not been discussed previously, such as the western side of the active region (the general 
area from x > 200 and y > 200). The phase speeds in the 3.5 mHz maps are of the order of several 
hundred km s _1 , which is to be expected for standing waves or propagating fast-mode waves in the 
corona. 

The 5.5 mHz maps (Figure ITOb are again much more sparsely populated, with only a small 
subset of locations showing significant oscillatory power. This is perhaps not surprising since 
previous analyses have shown that this particular flare led to transverse oscillations with periods of 
order five minutes. However, the fact that we do find power at 5.5 mHz (sithree minute periods), 
and that it occurs at mostly the same locations as in the 3.5 mHz maps may suggest that the 
transverse oscillations are not dominated by a single peak in the Fourier spectrum, but rather have 
a range of periods. 

The 1.5 mHz maps (Figure [8]) are very intriguing with most of the dimming regions showing 
significant oscillatory power that propagates at phase speeds of order 100-500 km s _1 . However, in 
this case it is more difficult to clearly distinguish whether the resulting phase speeds and locations 
are uniquely the sites of an oscillatory disturbance propagating, or whether the sudden coronal 
dimming triggered by the flare is partially responsible for the results. However, visual inspection 
of the Fourier-filtered time series suggests that the coronal dimming is indeed accompanied by 
transient and quasi-periodic "wiggling" of at least some of the coronal structures in these regions. 
The observed range of phase speeds may then well be the propagation speed of the fast mode waves 
that are associated with such a drastic reorganization of the corona. 

In the above we have focused on the locations of significant oscillatory power and the associ- 
ated phase speeds of the observed oscillatory disturbances. Our algorithm also naturally provides 
the direction in which the propagation occurs, as well as the coherence length and width of the 
islands of high coherence. The observed propagation direction usually is aligned with the domi- 
nant intensity structure (coronal loop or fan) that carries the disturbance. The coherence width and 
length are measures of how large the coherence island is, with the width by definition shorter than 
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the length, and the difference between both larger as the island of high coherence is more oblong 
in shape. Typical coherence lengths are of order three - six Mm in dataset I, with longer coherence 
lengths for lower frequencies. For dataset II coherence lengths are of similar size, except for the 
low frequency results which have coherence lengths that are up to 10 - 35 Mm in the regions where 
the coronal dimming dominates. This is not surprising, since the dimming impacts a large area of 
the active region. 



3.2. Reducing False Positives 

As we mentioned in the above, our algorithm is so sensitive that it finds a significant number 
of false positives - locations where instrumental or cosmic ray artifacts lead to spurious coherent 
signals. In addition, many of the locations have values for the phase speed that are based on poor 
fits between the distance along the propagation direction and the travel time, which can lead to 
significant errors on the parameters that are calculated. While the value of the maps we calculate 
is clear, the required human interpretation of what constitutes a false positive would render an 
automated classification (which is crucial for SDO/AIA) difficult. Fortunately we have found a 
method (which we describe in the following) to significantly reduce the number of false positives. 

We take several steps to prune the results down to locations where the observed signal is most 
likely caused by a real coronal signal, as opposed to an instrumental or detection artifact. We start 
with a map that contains locations that have enough neighbors (as defined in Sect. [B with high 
coherence. As a first step to winnow down the number of false positives, we use the error on 
the determination of the phase speed (from the least-squares fit of distance versus travel-times). 
We reject all locations for which the relative error on the phase speed is larger than 0.3 {i.e., the 
error divided by the phase speed). This removes from our map a large number of locations where 
the travel times are not well correlated with distance along the island of high coherence. Such 
a poor correlation leads to a poor linear fit and a large error on the phase speed. To reduce the 
large number of locations caused by instrumental artifacts we then reject all locations in which 
the average brightness is below a certain threshold (1 DN s _1 ). This removes a large number of 
artifacts that are caused by read-out noise (which can dominate coherent signals in the dark areas 
of the field of view), as well as artifacts from JPEG compression around cosmic rays. While the 
despiker usually removes the cosmic rays, the compression artifacts around them can sometimes 
remain in the data. Most of these artifacts are removed as a result of the brightness thresholding. 

As a next step we use the map of locations and label each contiguous region in the map using 
the IDL routine label_region.pro. We then reject all regions that contain fewer than, say, eigth 
pixels (superpixels of 1 x 1 arcsecond 2 in our case). This removes most of the moss regions and 
many of the instrumental artifacts that usually show coherence over only a few pixels. Note that 
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this step bases the rejection on the number of contiguous pixels in our map that have a significant 
number of neighbors that are highly coherent at the frequency studied. This a different from the 
step in Sect. \T\ where we rejected a location from the map if it had fewer than six neighboring 
pixels in the island of high coherence. While these two measures are related, they are not identical. 
The end result of these three steps is usually a cleaner map of locations with significant oscillatory 
power that are coronal in nature (as opposed to instrumental). We now have a map of regions (as 
opposed to pixels), in which adjacent pixels have been clustered together into one region. We then 
calculate the average phase speed, the average error on the phase speed, the average coherence 
length and width, and the average angle for each of these regions. The resulting figures show 
maps of oscillatory "events" with their averaged properties (see Figure [TT] and [14] for dataset 1 
and 2 respectively). The pruning method reduces the number of "pixels" from 850, 409, and 129 
respectively for 1.5, 3.5, and 5.5 mHz and combines them to 46, 14, and 6 regions in dataset I. 
Dataset II contains many more pixels with significant oscillatory power, but the pruning method 
reduces the number of pixels from 950, 1103, and 1150 respectively down to 91, 94, and 55 
events or regions. Such a significant reduction of information in an automated fashion renders the 
combination of our wavetracking and pruning algorithms into an interesting option to automatically 
mine the SDO/AIA data. The rejection of false positives is usually good enough so that the standard 
deviation of properties in each region is quite small (and is not particularly useful as a further 
method of rejecting artifacts). 

The rejection of false positives is generally a complicated task that is quite dependent on the 
quality of the data (exposure times, number of cosmic ray hits, efficiency of cosmic ray removal, 
etc.). As a result it usually requires some fine tuning of parameters. Fortunately the quality of the 
SDO/AIA data is expected to be higher and of a more constant quality than that of the TRACE data 
used here and so we will leave much of the "fine tuning" of the pruning thresholds, etc. until we 
have AIA data. This should greatly facilitate full and reliable automation of our data processing 
pipeline on the SDO/AIA data. 



4. Discussion 

It is clear from our results that our method reliably finds locations of significant oscillatory 
power in the corona. The two-step approach involves wavetracking and pruning out false-positives 
that can reliably and automatically detect events that are worthy of further detailed analysis. These 
events can be characterized by the dominant frequency, the average phase speed of the observed 
periodic disturbance, the average angle of propagation and the average coherence width and length 
(on a pixel by pixel basis), as well as the size of the cluster of coherent pixels. 

We have demonstrated that the phase speed of the propagating disturbance can be reliably used 
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to distinguish the various wave modes that are expected to appear in the corona, with low values 
(< 150 km s _1 ) of the phase speed associated with propagating slow-mode waves, and higher 
values (> 150 km s _1 ) associated with either propagating fast-mode waves or standing transverse 
oscillations. We note however that an accurate determination of the phase speed with the available 
data (and its relatively slow, and non-fixed cadence) is often difficult, so that the determined phase 
speeds are only indicative of the real value of the phase speed. Fixed high cadence, high S/N data 
such as that expected from SDO/AIA will alleviate much of this problem. In addition, we have 
experimented with the threshold of minimum coherence that defines the island of high coherence 
and find that higher thresholds generally of course reduce the size of the coherent islands, but can 
also lead to better determination of the phase speed. 

An additional problem that impacts the determination of the phase speed is the fact that we 
basically assume the island of coherence to be a roughly linear feature, such as a coronal loop. 
However, in certain cases, especially in fans associated with sunspots, or the large coronal dimming 
regions such an assumption may not be the most representative of how the propagation actually 
occurs. In such cases the source of the wave is relatively compact and the waves actually propagate 
outward in a more spherical manner, and not in a linear fashion, i.e., along a straight line which is 
what our method assumes. That means that our definition of distance travelled (which we define as 
the distance along the long axis of the island of high coherence) does not represent a good measure 
of distance along the propagation direction. This can clearly lead to a less accurate determination 
of the phase speed, and may be why some of the measured phase speeds of the slow-mode waves in 
the sunspot fan are on the low side (50 km s _1 ). This problem also directly affects the determination 
of the other parameters, especially the propagation angle and to a lesser extent the coherence width 
and length. Despite these issues, the phase speeds are different enough to distinguish propagating 
slow-mode waves from transverse oscillations. Our algorithm routinely finds phase speeds of order 
several hundred km s -1 for the latter, which is to be expected for partially standing waves. We note 
that previous analyses of these transverse oscillations have indicated that there is often a significant 
propagating component as well as a standing component to the oscillations. This may be one of 
the reasons why the phase speeds we obtain are not infinite. Other reasons, of course, are related to 
the inherent uncertainties in determining travel times and travelled distances, as described above. 

In its current form, our algorithm shows great potential for automated processing of SDO/AIA 
data. It is crucial for the SDO/AIA data that any kind of automated flagging of locations with sig- 
nificant oscillatory power contains so few false positives that the resulting list of events does not 
overwhelm the so-called "Heliophysics Knowledge Base" ( jhttp : / /www . lmsal . com/helio-inf ormati 
and at the same time significantly reduces the amount of information the end-user needs to sift 
through before finding the subset of data that contains the oscillations of interest. As we have 
shown in this work, our algorithm does exactly that. It provides of order a few to a few dozen 
events per hour of EUV data and allows the users to immediately home in on the regions of in- 
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terest that deserve more detailed study. Our algorithm takes one - two minutes (depending on the 
amount of oscillations present - we note that this does not include the time taken to derotate and 
coalign the image cubes) for a 256 x 256 x 1 10 datacube on a Mac Pro computer, without any effort 
to optimizing the FFT and correlation calculations for higher speeds. AIA data will be 2048x2048 
at this spatial resolution (i.e., two x two rebinning), so that a single processor would require of 
order one - two hours to process one hours worth of data in a single passband. This means that 
for all coronal AIA data, a pipeline containing of order twenty CPUs would be enough to keep 
pace with the data flow. The removal of false positives will require some finetuning of parameters 
and pruning approaches (e.g., using the cross-spectral power instead of the intensity to reject false 
positives) that will be investigated in the lead-up to the launch of SDO, and since much will depend 
on the quality of the data, during the commissioning phase. 



5. Summary 



Our algorithm was applied to two datasets that contain a diverse range of oscillation or wave- 
like features. We were able to identify moss oscillations mostly at 1.5-3.5 mHz, propagating slow- 
mode waves at 3.5-5.5 mHz, partially standing transverse oscillations at 1.5-5.5 mHz, as well as a 
large number of previously unreported apparently propagating waves at 1 .5 mHz in coronal loops 
associated with plage (dataset I). Our algorithm has the potential to perform automated studies of 
cross-wavelength cor relations of oscillatory power in various EUV wavelengths (171 and 195A 
for example, see, e.g. jKing et a/.ll2003l : iRobbrecht et a/.ll200l|) . This will be the subject of a future 
paper. Our algorithm can also be run on a large number of datasets to help determine how common 
propagating slow-mode oscillations are in coronal loops associated with sunspots and plage. Espe- 
cially the latter have been the subject o f speculation since they involve the leakage of a potentially 



sign ificant amount of p-mode power (|De Pontieu. Erdelvi. and De Moortel 



2006Q which could impact the amplitude of the p-mode spectrum (|de Moortel and Rosnenl2007T) . 



2005 



Jefferies et al. 



The presence of significant oscillatory power in the 1.5 mHz passband in non-flaring ac- 
tive region coronal loops provides us with some interesting scienti fic return. Previously some 



oscillatory power at ten minutes had be en detected in polar plumes (Deforest and Gurmanlll998 



Ofman. Nakariakov. and Deforest! 1 1999|) . but to our knowledge such power had not been reported 
in coronal loops in active regions. The signal is quite weak in these quiescent active-region loops, 
but most likely real. The fact that not all coronal loops show this propagating signal strongly 
suggests that slow co-alignment drifts are not the cause. We have tested this hypothesis by per- 
forming the same kind of calculations on dataset I using a variety of methods to remove the slow 
pointing drift that is intro duced into TRACE datasets from thermal flexing of the telescope tube 
(lAschwanden et a/.ll2000|) . The measured propagation speeds of these 1.5 mHz waves are of or- 
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der 50 km s -1 , similar to that of the well-known 3.5 mHz propagating slow-mode oscillations. 
The presence of these 1.5 mHz waves is intriguing in part because their source in the lower atmo- 
sphere may well be different from that of the more well-known 3.5 mHz waves. The leakage of 
p-modes through the lower atmosphere has been suggested as a source for the prevalence of five 



minu t e period waves in t he corona by a variety of authors (IDe Pontieu. Erdelyi. and De Moortel 



20051 : Ijefferies et q/J|2006|) . However, p-modes lack significant power at 1 .5 mHz to drive what we 



see. This suggests that a different source may be wo rking to produce t hese waves with ten minute 
periods in active region coronal loops. Recent work (IStraus et al. has seen the identification 

of horizontally propagating internal gravity waves at what is tho ught to be the b oundary between 
the photosphere and chromosphere. Long-standing theory (e.g.. ILi ghthilll 1 1 9 6 7|) has these waves 
coupling strongly to those that propagate along magnetic field lines into the outer atmosphere. A 
detailed statistical study of these 1 .5 mHz waves may well reveal more details about what the real 
source of these disturbances is and whether these acoustic gravity waves couple into the outer at- 
mosphere. We should note that at least in this one dataset the 1.5 mHz waves are absent from the 
sunspot fan, and appear only in a subset of plage-associated loops. We should note that it is also 
possible that these 1.5 mHz waves are actually changes in the coronal loop emission itself as a 
result of changes in heating or cooling of the loops. Detailed visual inspection of the unfiltered and 
filtered movies does not allow for an unequivocal identification of the dominant process underlying 
the signal our algorithm finds. We intend to investigate these 1.5 mHz signals further in follow-up 
work by studying a wide range of different TRACE datasets. Another new scientific finding of our 
numerical experiments is the presence of a significant amount of low-frequency oscillatory power 
in the wake of a coronal dimming. More detailed studies will be necessary to determine why these 
low frequencies are so prevalent in the dimming regions and whether this can be used to determine 
physical parameters of the dimming region itself. 



6. Conclusion 

We have considered the problem of automatically (and robustly) isolating and extracting in- 
formation about waves and oscillations observed in EUV image sequences of the solar corona with 
a view to near real-time application to data from the Atmospheric Imaging Array (AIA) on the 
Solar Dynamics Observatory (SDO). We have found that a simple coherence / travel-time based 
approach detects and provides a wealth of information on transverse and longitudinal wave phe- 
nomena in the test sequences provided by the Transition Region and Coronal Explorer (TRACE). 
The results of the search can be robustly "pruned" (based on diagnostic errors) to minimize false- 
detections such that the remainder provides reliable measurements of waves in the solar corona, 
with the calculated propagation speed allowing automated distinction between various wave modes 
and can be automatically applied to the enormous flow of data (^lTb day -1 ) that will be provided 
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by SDO/AIA. In addition to demonstrating the applicability of this approach we have found signa- 
tures of some interesting low frequency (1.5 mHz) oscillatory phenomena in the TRACE test data 
that motivate further study. 

SWM is supported by NSF ATM-0541567, NASA NNG06GC89G; BDP by NASA grants 
NAS5-38099 (TRACE), NNM07AA01C (Hinode) and NNG06GG79G. SWM and BDP are jointly 
supported by NASA grants NNX08AL22G and NNX08AH45G. 
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Fig. 1. — Examples of the coherence, weighted cross-power, phase difference and phase travel- 
time for for the 1 July 1998 (panels A-D) and 14 July 1998 (panels E - H) datasets for 5.5 mHz 
and 1 .5 mHz Fourier filter samples respectively. In each case the reference pixel is at 0,0, the solid 
black contour outlines regions of highly coherent signal and the fitted angle (solid white lines in 
panels C and G) to central coherence "island" is shown as a solid white line in panels C (14.6°) 
and G (-30.4°). 
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Fig. 2. — Examples of phase speed determination for the cases shown in Figure [Q In each case 
the green triangles mark the positions from the reference pixel and measured phase travel-time at 
that pixel. The gradient of the least-squares linear fit (red dot-dashed line) yields the phase-speed 
of the propagating signal. 



-18- 




Fig. 3. — Example pixel-to-pixel travel-time analysis in the temporal domain. In panels A and B 
we show the TRACE 171 A lightcurve (black solid line) and its 5.5 mHz filtered counterpart (blue 
solid line) between the reference pixel [181, 97] and pixel [179, 97] for the 1 July 1998 dataset (cf. 
Figured)). In panel C we show the derived cross-correlation function (black solid line), envelope 
maxima (red triangles) and fit to the central peak (green solid line). The vertical dashed lines 
in panel C show the positions of the estimated group travel-time (red; -4.17 minutes) and phase 
travel-time (green; -0.527 minutes, i.e., consistent with the values in Figures[T]and propagating left 
to right which is consistent with the movie of Figure respectively. 
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Fig. 4. — Illustrations of the two TRACE 171 A timeseries used in the presented analysis; 1 July 
1998 (left) and 14 July 1998 (right). In each panel we show (clockwise from the top left) snapshots 
of the TRACE 171 A intensity, 1.5 mHz, 5.5 mHz, and 3.5 mHz filtered timeseries respectively. The 
boxes marked on the panels correspond to the regions used to construct Figure[T]and the subsequent 
analyses. See the online edition of the Journal for animations of these figures. 
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Fig. 5. — Results of the coronal wave detection algorithm for the 1.5 mHz filtered timeseries for 
the 1 July 1998 dataset. We show the average TRACE 171 A intensity image, weighted signal 
coherence, length, width, angle, phase speed, and relative error of the phase speed computed from 
the large coherence island. 
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Fig. 6. — Results of the coronal wave detection algorithm for the 3.5 mHz filtered timeseries for 
the 1 July 1998 dataset. We show the average TRACE 171 A intensity image, weighted signal 
coherence, length, width, angle, phase speed, and relative error of the phase speed computed from 
the large coherence island. 
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Fig. 7. — Results of the coronal wave detection algorithm for the 5.5 mHz filtered timeseries for 
the 1 July 1998 dataset. We show the average TRACE 171 A intensity image, weighted signal 
coherence, length, width, angle, phase speed, and relative error of the phase speed computed from 
the large coherence island. 
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Fig. 8. — Results of the coronal wave detection algorithm for the 1.5 mHz filtered timeseries for 
the 1 July 1998 dataset. We show the average TRACE 171 A intensity image, weighted signal 
coherence, length, width, angle, phase speed, and relative error of the phase speed computed from 
the large coherence island. 
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Fig. 9. — Results of the coronal wave detection algorithm for the 3.5 mHz filtered timeseries for 
the 14 July 1998 dataset. We show the average TRACE 171 A intensity image, weighted signal 
coherence, length, width, angle, phase speed, and relative error of the phase speed computed from 
the large coherence island. 




Fig. 10. — Results of the coronal wave detection algorithm for the 5.5 mHz filtered timeseries for 
the 14 July 1998 dataset. We show the average TRACE 171 A intensity image, weighted signal 
coherence, length, width, angle, phase speed, and relative error of the phase speed computed from 
the large coherence island. 
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Fig. 11. — Final results of the coronal wave detection algorithm for the 1.5 mHz filtered timeseries 
for the 14 July 1998 dataset that have been "pruned" to reduce the appearance of false positive 
detections using the technique discussed in Sect. 13.21 We show the average TRACE 171 A intensity 
image and the region averaged signal coherence, length, angle and phase speed, and the errors in 
the latter. 
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Fig. 12. — Final results of the coronal wave detection algorithm for the 3.5 mHz filtered timeseries 
for the 1 July 1998 dataset that have been "pruned" to reduce the appearance of false positive 
detections using the technique discussed in Sect. 13.21 We show the average TRACE 171 A intensity 
image and the region averaged signal coherence, length, angle and phase speed, and the errors in 
the latter. 
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Fig. 13. — Final results of the coronal wave detection algorithm for the 5.5 mHz filtered timeseries 
for the 1 July 1998 dataset that have been "pruned" to reduce the appearance of false positive 
detections using the technique discussed in Sect. 13.21 We show the average TRACE 171 A intensity 
image and the region averaged signal coherence, length, angle and phase speed, and the errors in 
the latter. 
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Fig. 14. — Final results of the coronal wave detection algorithm for the 1.5 mHz filtered timeseries 
for the 14 July 1998 dataset that have been "pruned" to reduce the appearance of false positive 
detections using the technique discussed in Sect. 13.21 We show the average TRACE 171 A intensity 
image and the region averaged signal coherence, length, angle and phase speed, and the errors in 
the latter. 
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Fig. 15. — Final results of the coronal wave detection algorithm for the 3.5 mHz filtered timeseries 
for the 14 July 1998 dataset that have been "pruned" to reduce the appearance of false positive 
detections using the technique discussed in Sect. 13.21 We show the average TRACE 171 A intensity 
image and the region averaged signal coherence, length, angle and phase speed, and the errors in 
the latter. 
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Fig. 16. — Final results of the coronal wave detection algorithm for the 5.5 mHz filtered timeseries 
for the 14 July 1998 dataset that have been "pruned" to reduce the appearance of false positive 
detections using the technique discussed in Sect. 13.21 We show the average TRACE 171 A intensity 
image and the region averaged signal coherence, length, angle and phase speed, and the errors in 
the latter. 
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